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Abstract 


This  paper  describes  an  algorithm  for  reducing  a real 
matrix  A to  block  diagonal  form  by  a real  similarity  trans- 
formation. The  columns  of  the  transformation  corresponding 
to  a block  span  a reducing  subspace  of  A,  and  the  block  is 
the  representation  of  A in  that  subspace  with  respect  to 
the  basis.  The  algorithm  attempts  to  control  the  condition 
of  the  transformation  matrices,  so  that  the  reducing  sub- 
spaces are  well  conditioned  and  the  basis  vectors  are  numeri- 
cally independent. 
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AN  ALGORITHM  FOR  COMPUTING  REDUCING  SIJBSPACES 
BY  BLOCK  DIAGONALIZATION 


Connice  A.  Bavely 
G.  W.  Stewart 


1.  Introduction 

The  purpose  of  this  report  is  to  describe  an  algorithm  for  reducing 
a real  matrix  A of  order  n to  a block  diagonal  form  by  a real  simi- 
larity transformation.  Specifically,  the  algorithm  attempts  to  compute 
a real  nonsingular  matrix  X such  that  X ^AX  has  the  form 

(1.1)  X'^AX  = B = diag(B^,B2 B^)  , 

where  each  matrix  B^  is  square  of  order  n^. 

A deconposition  such  as  (1.1)  has  many  applications.  When  the  blocks 
B^  are  small,  powers  of  A can  be  economically  calculated  in  the  form 

a’^  = X diag(B^,B^,...,B^)X‘^  , 

and  this  fact  can  be  used  to  simplify  the  computation  of  functions  of  A 
defined  by  power  series  (e.g.  see  [7]).  If  X is  partitioned  in  the  form 


X = (Xj,X2,...,X3)  , 

where  each  X^  has  n^  columns,  then  the  columns  of  X^  form  a basis  for 
a reducing  subspace  of  A,  and  B^  is  the  representation  of  A with  respect 
to  that  basis.  The  associated  spectral  projector  is  given  by  where 
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is  formed  from  the  corresponding  rows  of  X ^ (for  definitions 
and  applications  see  [4]). 

There  are  theoretical  and  practical  limitations  on  how  small  the 
blocks  in  (1.1)  can  be.  Theoretically,  they  can  be  no  smaller  than  the 
blocks  in  the  Jordan  canonical  form  of  A.  Practically,  they  may  have 

to  be  larger.  The  numerical  problems  associated  with  decompositions  such 
as  (1.1)  have  been  examined  in  detail  in  [3].  Here  we  give  only  a brief 
sutimary. 

The  principal  difficulty  is  that  the  Jordan  form  of  a matrix  need  not 
be  numerically  well  determined;  very  small  perturbations  in  the  matrix 
may  cause  blocks  to  split  or  coalesce.  Any  atten^Jt  to  separate  two  such 
"nearby"  blocks  will  result  in  a transformation  matrix  X whose  columns 
are  nearly  linearly  dependent,  or  equivalently  X will  be  ill-conditioned 
in  the  sense  that  the  product  ||X||||X  ^||  is  large  (here  H*!!  denotes  a 
suitable  matrix  norm).  In  this  case,  it  will  be  iir^ossible  to  form  X ^ 
or  solve  linear  systems  involving  X accurately  [10,14].  The  phenomenon 
is  loosely  associated  with  close  eigenvalues;  but  there  are  matrices  with 
equal  eigenvalues,  e.g.  symmetric  matrices,  that  can  be  split  completely 
into  1x1  blocks,  and  there  are  matrices  with  well  separated  eigenvalues 
that  cannot  be  split  except  by  very  ill-conditioned  transformations. 

Our  algorithm  attempts  to  avoid  these  difficulties  by  working  only 
with  well-conditioned  transformations.  If  a group  of  eigenvalues  cannot 
be  split  off  into  a block  by  a transformation  whose  condition  observes  a 


tolerance  provided  by  the  user,  the  block  is  enlarged  until  a well- 
conditioned  reducing  transformation  can  be  found.  In  principle  this  does 
not  insure  that  the  final  transformation  will  be  well -conditioned,  since 
it  is  formed  as  the  product  of  a number  of  reducing  transformations; 
however,  we  have  found  that  when  a matrix  possesses  a well-conditioned 
decOTiposition  of  the  form  (1.1),  our  algorithm  generally  finds  it.  And 
the  exceptions  have  not  so  much  to  do  with  the  ill -conditioning  of  X as 
with  the  failure  of  the  algorithm  to  split  the  matrix  completely  owing 
to  the  comingling  of  degenerate  eigenvalues  with  well-conditioned  ones. 

A good  deal  of  work  has  been  done  on  the  numerically  stable  siit^lifi- 
cation  of  matrices  by  similarity  transformations  [5,6,8,12,13],  most  of 
which  has  been  sunmarized  in  [3].  For  the  most  part,  these  algorithms 
atten5)t  to  go  farther  than  ours  in  reducing  the  matrix,  however  at  con- 
siderable cost  in  complexity  and  computation.  The  virtues  of  the  algo- 
rithm proposed  here  are  its  sijiplicity  and  economy.  IVhen  it  is  required 
to  reduce  a matrix  beyond  what  is  done  by  our  algorithm,  the  other 
techniques  can  be  applied  to  the  blocks  produced  by  our  algorithm.  The 
algorithm  also  has  the  advantage  that  it  works  entirely  with  real  matrices 
by  the  device  of  groiping  pairs  of  coitplex  conjugate  eigenvalues  in  the 
same  block. 

In  the  next  section  of  this  paper  the  algorithm  is  described.  In 
Section  3 scrnie  numerical  exanples  are  given.  Programming  details  and  a 
listing  are  given  in  an  appendix. 
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2.  The  algorithm. 

The  first  part  of  the  algorithm  uses  orthogonal  transformations  to 
reduce  the  matrix  A to  quasi -triangular  form,  that  is  to  a block  iqjper- 
triangular  form  in  which  the  diagonal  blocks  are  of  order  at  most  two. 

The  blocks  of  order  one  contain  real  eigenvalues  of  A and  the  blocks 
of  order  two  contain  complex  conjugate  pairs  of  eigenvalues.  The  order- 
ing of  the  blocks  is  arbitrary,  and  the  order  can  be  changed  by  applying 
appropriate  orthogonal  transformations.  Since  this  reduction  of  A can 
be  effected  by  standard  techniques  [9,11] , we  may  assume  that  A is 
already  in  quasi-triangular  form. 

The  subsequent  block  diagonal ization  is  acconplished  as  follows. 

The  matrix  A is  partitioned  in  the  form 


A = 


f 


where  initially  A^^^^  is  1 x i or  2 x 2 depending  on  the  dimension  of 
the  leading  diagonal  block  of  A.  An  attenpt  is  then  made  to  find  a 
similarity  transformation  X such  that 


If  such  a transformation  can  be  found  and  if  it  is  not  too  ill-conditioned, 
the  reduction  proceeds  with  the  submatrix  A22*  If  not,  a suitable  1 x l 
or  2 X 2 diagonal  block  from  A22  is  located  and  moved  by  means  of  ortho- 
gonal transformations  to  the  leading  position  of  A22*  The  block  is  then 


adjoined  to  increasing  the  order  of  by  one  or  two,  as  is 

appropriate,  and  another  attempt  is  made  to  find  a reducing  matrix  X. 


The  iji^jlementation  of  such  an  algorithm  requires  the  answers  to  two  i 

qiiestions.  j 


1.  How  may  the  transformation  X be  ccmqjuted? 

2.  In  the  event  of  failure,  which  block  of  A22  is 


I 

f 


to  be  incorporated  into  j 

i 

We  shall  now  answer  these  questions.  | 

i 

We  seek  the  transformation  X in  the  form  i 

! 


(2.1) 


where  the  identity  matrices  are  of  the  same  orders  as  A^j  and  A22*  The 
inverse  of  X is  easily  seen  to  be 


(2.2) 

Hence 


X'^AX 


A^l  A^iP*PA22'^Aj^2 


and  the  problem  of  determining  X becomes  that  of  solving  the  equation 


(2.3) 


"^11^  ■ ^^22 


Because  A^^^  and  A22  are  quasi -triangular,  this  equation  can  be  solved 
by  a back- substitution  algorithm  of  Bartels  and  Stewart  [1] , provided  the 


eigenvalues  of  and  A22  are  disjoint. 

From  (2.1)  and  (2.2)  it  follows  that  X will  be  ill-conditioned 
whenever  P is  large.  As  each  element  of  P is  generated,  it  is  tested 
to  see  if  its  magnitude  exceeds  a bound  provided  by  the  user.  If  it  does, 
the  attenpt  to  confute  X is  abandoned  and  a new,  larger  block  A^j  is 
formed.  If  no  element  of  P exceeds  the  bound,  the  matrix  X is  ac- 
cepted and  the  matrix  A is  deflated  as  described  above.  The  transforma- 
tion X is  postmultiplied  into  a matrix  that  accumulates  all  the  trans- 
formations made  on  the  matrix  A. 

The  process  for  selecting  alxlorZxZ  block  of  A22  bo  incor- 
porate into  Aj^j^  goes  as  follows.  We  conpute  the  mean  of  those  eigen- 
values of  Aj^j^  having  nonnegative  imaginary  part.  A block  is  chosen 
from  A^2  whose  eigenvalue  with  nonnegative  imaginary  part  is  nearest 
this  mean.  This  block  is  moved,  as  described  above,  by  orthogonal  trans- 
formations to  the  leading  position  in  A22»  where  it  is  incorporated  into 
The  program  HQR3  [11],  which  can  be  used  to  obtain  the  initial  quasi- 
triangular  form,  has  a subroutine  which  will  compute  these  orthogonal 
transformation.  The  transformations  are  of  course  postmultiplied  into  the 
accumulating  matrix. 

We  sumnarize  our  algorithm  in  the  following  informal  code.  Further 
details  can  be  found  in  the  appendix  to  this  paper,  where  a FORTRAN  sub- 
routine in^jlementing  this  code  is  given.  The  code  takes  as  input  an 
array  A of  order  N containing  the  matrix  A and  an  array  X in  which 
the  transformations  are  accumulated.  In  addition  the  user  must  provide  a 


■i  in»*ii  ■ 1^1  tiiifi  .ilfc  ~ ' 


I 

i 

i 
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tolerance  to  bound  the  size  of  the  elements  of  the  deflating  transforma- 
tions. The  integers  Lll  and  L22  point  to  the  beginnings  of  the  current 
blocks  Ajj  and  A22  in  the  array  A.  The  informal  code  should  be  self- 
explanatory.  Comments  are  delineated  by  the  symbol  #. 


2 

3 

3.1 

3.2 

3.3 
3.3.1 
3.3. It. 1 
3. 3. It. 2 
3.3.1 
3.3.1e.l 

3.3.1e.2; 


3.3.1e.3: 

3.3.1e.4 

3.3.1 

3.3.2 

3.3.3 

3.3.4 

3.3.5 

3.3 

3.4 


reduce  A to  quasitriangular  form,  accumulating  the 
transformations  in  X [10]; 

Lll  = 1; 

loop  # until  the  matrix  is  diagonalized  # 
if  Lll  > N then  leave  3 fi; 

IJ2  = Lll; 

loop  # until  a block  has  been  deflated  # 

if  L22  = Lll  then  # use  the  first  1 x 1 or  2 x 2 block  # 

M = order  of  the  block  at  Lll; 

L22  = L22  + M; 

else  # augment  A, , with  alxlor  2 x 2 block  from  A,,  # 
confute  the  mean'^of  the  eigenvalues  of  A,  with 
nonnegative  imaginary  parts; 

find  the  MxM#M=lor2#  block  of  A22  whose  eigen- 
value with  nonnegative  imaginary  part  is^nearest  the 
mean; 

move  the  block  to  the  leading  position  of  A22  accumulating 
the  transformations  in  X; 

L22  = L22  + M # which  incorporates  the  block  in  #; 

fi; 

IT  L22  > N th^  leave  3.3  fi; 
attempt  to  split  off  A, , Tl  ] ; 
if  the  attempt  was  successful  then  leave  3.3  fi; 
restore  A,-; 
end  loop; 

if L22  £ N then  accumulate  the  deflating  transformation  in  X 

scale  columns  Lll  through  L22-1  of  X so  that  they  have 
2-norm  unity,  adjusting  A accordingly; 

Lll  = L22; 
end  loop; 


Several  comments  should  be  made  about  the  algorithm.  First,  it  uses 
only  real  arithmetic,  even  when  A has  complex  eigenvalues.  Second,  the 
algorithm  cannot  guarantee  that  the  final  transformation  is  well  condi- 
tioned, since  the  bound  on  the  elements  of  P restricts  the  condition  of 
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only  the  individual  transformations  craiprising  the  final  one.  Nonethe- 
less, we  have  found  little  tendency  toward  excessive  growth  in  the  condi- 
tion of  the  transformation.  TTiird,  no  att&npt  is  made  to  segregate  nearly 
equal  eigenvalues  initially  into  clusters;  whether  an  eigenvalue  splits 
off  or  not  depends  entirely  on  its  determining  a suitably  small  matrix 
P.  This  is  important,  since  it  means  that  the  algorithm  can  compute 
well  conditioned  eigenvectors  for  multiple  eigenvalues  (a  little  help  is 
required  from  rounding  error;  see  Section  3). 

The  strategy  for  determining  what  eigenvalues  to  add  to  the  current 
group  has  the  defect  that  it  can  mix  well-conditioned  and  ill-conditioned 
eigenvalues  that  are  nearly  equal,  thus  missing  the  possibility  of  a more 
conplete  reduction.  This  is  not  a serious  problem  when  the  blocks  are 
small.  However,  if  a finer  resolution  of  the  structure  of  the  matrix  is 
required,  the  techniques  referenced  in  Section  1 may  be  applied  to  the 
blocks  produced  by  our  algorithm.  In  fact  our  algorithm  can  be  regarded 
as  a preprocessing  step  to  reduce  the  problem  to  a size  where  it  can  be 
attacked  by  more  sophisticated,  but  more  expensive  methods. 

We  note  that  the  output  of  our  algorithm  may  depend  on  the  user 
svqpplied  tolerance  for  the  elements  of  P.  In  general  the  larger  the 
tolerance,  the  smaller  the  blocks  but  the  more  ill  conditioned  the  trans- 
formation. This  tradeoff  is  an  inevitable  consequence  of  the  poor  deter- 
mination of  the  structure  of  a matrix  in  the  presence  of  errors,  and  we 
know  of  no  algorithm  that  relieves  the  user  of  the  necessity  of  making  a 
decision  of  this  kind. 


2 

So  fax'  as  storage  is  concerned,  the  algorithm  requires  2n  loca- 
tions to  contain  the  matrices  A and  X and  a number  of  arrays  of  order 
n.  This  is  the  same  as  the  storages  required  by  algorithms  that  compute 
the  eigenvectors  of  a general  matrix. 

Excluding  the  initial  reduction  of  A to  quasitriangular  form,  the 

bulk  of  the  canputations  in  the  algorithm  occur  at  statement  3.3.3,  where 

an  attenqit  is  made  to  solve  the  equation  (2.3),  and  at  statement  3.4,  where 

the  transformation  is  accumulated.  The  multiplication  coimt  for  the  algo- 

2 2 

rithm  for  solving  (2.3)  is  of  order  (il  m + m Jl)/2,  where  K.  is  the  size 
of  Ajj^  and  m is  the  size  of  A22*  The  cost  of  accumulating  this 
transformation  in  X is  of  order  i,*m*n.  Thus,  at  one  extreme,  if  all  the 
eigenvalues  of  A are  real  and  they  all  can  be  deflated,  the  cost  in 
multiplications  will  be  of  order  n^/2,  which  compares  favorably  with  algo- 
rithms for  confuting  the  eigenvectors  of  A from  its  quasitriangular 
form.  On  the  other  hand,  if  the  eigenvalues  of  A are  real  and  none  of 
them  can  be  deflated,  the  algorithm  will  require  on  the  order  of  n^/12 
multiplications. 

4 

Although  an  n operation  count  is  disturbing,  there  are  several 
mitigating  factors.  First,  we  do  not  expect  the  algorithm  to  find  many 
applications  to  matrices  that  cannot  be  substantially  reduced  by  it, 
since  the  object  of  using  it  is  to  save  on  subsequent  calculations  with 
the  matrix.  Second,  the  count  assumes  that  the  algorithm  for  solving  (2.3) 
must  be  executed  to  completion  before  it  is  found  that  P is  unacceptably 
large.  This  is  not  likely;  in  fact  because  the  algorithm  [11]  for  reducing 
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A to  quasitriangular  form  arranges  the  eigenvalues  in  decreasing  order 
of  magnitude,  it  is  rather  unlikely.  Finally  the  order  constant  1/12  is 
is  modest;  for  n less  than  60  the  bound  is  less  than  5n^. 


3.  Numerical  results. 

In  this  section  we  summarize  the  results  of  some  numerical  experi- 
ments. Two  of  the  tests  were  performed  with  a class  of  test  matrices 
generated  as  follows.  Let  J be  a given  matrix  whose  structure  is  known 
(e.g.  J could  be  in  Jordan  canonical  form).  Let 


2 T 

Hi  = I - i ee^ 


vdiere  e = (1,1,...,!)^,  and 

H-  = I - - ff^  , 

2 n ’ 

where  f = (1,-1,. .. ,(-!)”  ) . Note  that  and  H2  are  symmetric  and 
orthogonal . Let 


S = diag(l,a,a^, . . . ,a^  ^) 


vdiere  a is  a given  parameter.  Then  we  take  A in  the  form 
(3.1)  A=  (H2SH^)J(H2SH^)'-'  ==  (H^SH^ J(H^S‘^H2)  . 


The  matrix  A,  vdiich  is  cheap  to  compute,  has  the  same  structure  as  J. 
The  transformation  H2SHJ  can  be  made  as  ill-conditioned  as  desired  by 
varying  a. 


J 
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In  describing  the  results  of  the  tests  we  report  two  nmbers.  The 
first  is  the  maximcn  p of  the  scaled  residuals 

1IAX.-XB.IL 

"i'  WI.IPtilL  ’ 

\diere  is  conposed  of  the  columns  of  X corresponding  to  the  block 
Second  we  report  ||X*^(L.  Together  these  numbers  give  an  idea  of  how 
stably  the  reduction  proceeded;  for  if 

R = AX  - XB 

then,  with  E = RX  we  have  that 

(A+E)X  - XB  = 0 ; 

that  is  B is  exactly  similar  to  A + E.  The  relative  error  that  this  per- 
turbation represents  in  A is 


Since  ||X||  1,  the  relative  error  will  be  of  the  order  p(|X  ^|( . 

The  first  test  case  illustrates  the  ability  of  the  algorithm  to  split 
apart  nearly  equal  eigenvalues  with  independent  eigenvectors.  We  took 

J = diagj^l,l-e,l+e,  Q J ^ ,.3,.4,.5,.6,.7]  . 

The  algorithm  was  applied  for  various  values  of  e,  a,  and  RMAX,  the  bound 
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on  the  size  of  the  deflating  transformations.  The  results  are  sumnarized 
in  Table  1,  in  which  the  eigenvalues  of  A are  numbered  as  follows: 

123456789  10 

1+i  1-i  1+e  1 1-e  .7  .6  .5  .4  .3 

Complex  eigenvalues  are  indicated  by  a circumflex. 

Provided  RMAX  was  large  enough  to  allow  a sufficiently  ill  conditioned 
transformation,  all  the  cases  were  split  completely.  The  condition  of  the 
transformation  was  in  no  case  much  greater  than  the  condition  of  the  scramb- 
ling transformation  in  (3.1).  It  is  of  interest  to  note  that  the  algorithm 
was  successful  when  e = 0,  that  is  when  A has  three  equal  eigenvalues. 
Mathematically,  the  algorithm  for  solving  (2.3)  breaks  down  when  A^^^ 
and  A22  have  common  eigenvalues;  however,  the  experiments  indicate  that 
if  the  equal  eigenvalues  have  independent  vectors,  rounding  error  will 
perturb  them  enough  for  the  algorithm  to  work. 

The  second  example  shows  the  failure  of  the  algorithm's  strategy  for 
selecting  the  next  eigenvalue  to  be  adjoined  to  A^^.  Here  J has  the 
form 


r / ^ 

1 

0 \ 

J = diag 

1,  0 

1 

1 j , . 3, .4, .5, .6, . 7, .8 

L \0 

0 

1 / J 

and  a was  taken  to  be  one.  Of  the  four  eigenvalues  at  unity,  one  is  per- 
fectly conditioned,  and  the  other  three,  which  belong  to  a single  Jordan 
block,  are  very  ill-conditioned.  To  six  figures  the  computed  eigenvalues 
were 
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1.00073  + 0.0013031 
1.00073  - 0.0013031 

1.00000 
0.99853 
.80000 

ITie  pair  of  conplex  eigenvalues  could  not  be  deflated,  since  they  were 
coupled  to  the  third  member  of  the  block.  But  this  member  would  not  be 
adjoined  without  first  adjoining  the  well  conditioned  eigenvalue  at  unity. 
Consequently,  the  algorithm  produced  a single  block  of  order  four,  rather 
than  two  blocks  of  orders  one  and  three.  This  block  of  order  four  could 
be  reduced  further  by  the  more  sophisticated  techniques  described  in  the 
references . 

The  third  test  case  is  the  Frank  matrix  of  order  12  which  has  appeared 
frequently  in  the  literature  [2,3].  The  smaller  eigenvalues  of  this  matrix 
are  extranely  ill  conditioned.  The  results  of  test  runs  on  this  matrix 
are  summarized  below. 


RMAX 

BLOCK  STRUCTURE 

P 

llx'^IL 

10 

1,2,3,4,5,(6,7,8,9,10,11,12) 

2.9  X 10'^ 

93.6 

50 

1,2,3,4,5,6,(7,8,9,10,11,12) 

2.9  X 10’^ 

2920.2 

100 

1,2,3,4,5,6,(7,8,9,10,11,12) 

2.9  X 10’^ 

2920.2 

1000 

1,2,3,4,5,6,(7,8,9,10,11,12) 

2.9  X 10"^ 

2920.2 

The  algorithm  performed  much  as  expected,  separating  the  larger  eigenvalues 
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and  groi5)ing  the  smaller  eigenvalues  together.  This  grouping  is  consistent 
with  the  precision  of  the  ccanputation. 

The  final  test  case  is  included  because  the  failure  of  our  algorithm 
to  decon^jose  it  reveals  the  shaky  foundations  of  a fairly  common  numerical 
practice.  Specifically  we  generated  the  conpanion  matrix  of  the  polynomial 
given  in  [14,  p.  74].  Since  the  zeros  of  this  polynomial  are  not  very  ill 
conditioned,  we  were  surprised  when  the  algorithm  failed  to  split  off  so 
much  as  one.  Some  further  conqjutations  revealed  that  the  matrix  of  left 
eigenvectors  (the  inverse  Vandermode  of  the  zeros)  had  rows  of  order 

f\  ft 

10°  - 10°.  Hiis  explains  the  failure  to  reduce  the  matrix.  More  inportant, 
though,  it  shows  that  the  eigenvalues  of  the  conpanion  matrix  are  much 
more  ill  conditioned  than  the  zeros  of  the  polynomial  and  suggests  that  the 
practice  of  using  eigenvalue  routines  to  find  zeros  of  polynomials  can 
result  in  an  unnecessary  loss  in  accuracy. 


a 

e 

RMAX 

BLOCK  STRUCTURE 

P 

oo 

1.0 

10’^ 

10. 

(1,2), 3,4, 5, 6, 7, 8, 9. 10 

2.1 

X 

10'^ 

7.7 

10'^ 

10. 

(1,2), 3, 4, 5, 6, 7, 8, 9, 10 

1.4 

X 

o 

1 

6.7 

0.0 

10. 

(1,2), 3, 4, 5.6, 7, 8, 9, 10 

1.1 

X 

10'^ 

8.0 

10.0 

10‘^ 

10. 

(1,2), 3, 4, 5, 6, 7, 8, 9,10 

1.5 

X 

10‘^ 

25.6 

10’^ 

10. 

(i, 2), 3, 4, 5, 6, 7, 8, 9, 10 

1.8 

X 

10'^. 

25.4 

0.0 

10. 

(1.2), 3, 4,5, 6,7, 8, 9, 10 

1.6 

X 

10'" 

24.2 

100.0 

o 

1 

10. 

(1,2), 3, 4, 5, 6, 7, 8, 9, 10 

1.1 

X 

10'^ 

187.9 

10'^ 

10. 

(1,2), 3, 4, 5, 6, 7, 8, 9, 10 

1.3 

X 

10'^ 

196.6 

0.0 

10. 

(1,2), 3, 4, 5. 6, 7, 8, 9, 10 

1.2 

X 

10'^ 

157.9 

1000.0 

10’^ 

10. 

(i, 2, 4,3, 5), (6, 7, 8, 9), 10 

1.1 

X 

10"^ 

326.6 

100. 

(1,2), 3, 4, 5, 6,7, 8, 9, 10 

1.0 

X 

10-« 

1537.9 

10'^ 

10. 

(1,2,4, 3, 5), (6, 7, 8, 9), 10 

7.8 

X 

10-s 

333.8 

100. 

(1,2), 3, 4, 5, 6, 7, 8, 9, 10 

7.8 

X 

10-8 

1559.1 

0.0 

10. 

(1,2,3), 4, 5, (6, 7, 8, 9), 10 

8.2 

X 

10-8 

394.1 

100. 

(1,2), 3, 4, 5, 6, 7, 8, 9, 10 

8.2 

X 

10-8 

1355.7 

Table  4.1 
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Appendix 

PROGRAMING  DETAILS  AND  PROGRAM  LISTING 

Al.  Usage. 

The  calling  sequence  for  BDIAG  is 

CALL  BDIAG(A,LDA,N,EPSHQR,RMAX,ER,EI,TVPE,BS,X,LDX,FAIL) . 

The  parameters  in  the  calling  sequence  are  (starred  parameters  are 


altered  by  the  subroutine) 


*ACLDA,N) 

* 

ar  array  that  initially  contains  the  N x N matrix 
to  be  reduced.  On  return  A contains  the  reduced 
block  diagonal  matrix. 

. LDA 

the  leading  dimension  of  A. 

N 

the  order  of  the  matrices  A and  X. 

EPSHQR 

a real  variable  containing  a convergence  criterion  for 
the  siibroutines  HQR3  and  EXCHNG  [11]. 

RMAX 

1 

a real  variable  containing  a bound  on  the  absolute 
values  of  the  elements  in  the  reducing  matrices. 

, *ER(N) 

1 

a real  array  containing  the  real  parts  of  the  eigen- 
values of  A. 

1 *EI  (N) 

a real  array  containing  the  imaginary  parts  of  the 
eigenvalues . 

*TYPE(N) 

1 

an  integer  array  whose  i-th  entry  is 

0 if  the  i-th  eigenvalue  is  real 

1 if  the  i-th  eigenvalue  is  complex  with  positive 

imaginary  part 

2 if  the  i-th  eigenvalue  is  complex  with  negative 

imaginary  part 

-1  if  the  i-th  eigenvalue  could  not  be  confuted 

1 *BS(N) 

a singly  subscripted  array  that  contains  infoimation 
on  the  block  structure  of  the  matrix  returned  by  the 
program.  If  there  is  a block  of  order  K at  A(L,L) , 
then  BS(L) ,BS(L+1) , . . .BS(L+K-1)  contain  the  integers 
K,- (K-1) , . . . ,-l.  Thus  a positive  entry  of  K indi- 
cates the  start  of  a block  of  order  K. 

7- 

* 

I'i 

i t 
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*X(LDX,N)  an  array  into  which  the  reducing  transformations  are 
accumulated. 

LDX  the  leading  dimension  of  X. 

FAIL  a logical  variable  which  is  false  on  a normal  return 

and  is  true  on  return  if  an  error  has  occurred. 

The  program  requires  the  programs  ORTHES  [9],  ORTRAN  [9],  HQR3  [11], 
EXCHNG  [11],  SHRSLV  [1],  DAD  and  SPLIT  [11]. 

A suitable  choice  for  the  parameter  EP9IQR  is  the  rounding  unit  of 
the  conputer  on  which  the  program  is  being  run;  i.e.  if  one  is  working 
in  a precision  of  about  t decimal  figures  EPSHQR  should  be  of  order  10 
A2.  Programming  details. 

In  this  section  we  shall  describe  sane  of  the  details  of  the  imple- 
mentation of  the  algorithm.  Throughout  this  section  we  refer  to  the  outline 
of  the  algorithm  in  Section  2. 

Statement  1.  This  is  accomplished  by  the  EISPACK  routines  ORTHES  and 
ORTRAN  [ 9]  and  the  QR  routine  HQR3  [11]. 

Statement  2.  Lll  points  to  the  leading  position  of  the  current 
block  All,  which  is  of  order  DAll. 

Statement  3.  This  is  the  main  loop  of  the  program.  It  ends  vdien 
Lll  > N,  indicating  that  there  are  no  further  blocks  to  deflate. 

Statonent  3.2.  L22  points  to  the  leading  position  of  the  current 

block  A22,  which  is  of  order  DA22. 

Statement  3.3.  In  this  loop  Aj^^^  is  repeatedly  augmented  until  it 
can  be  deflated  or  until  A22  is  void. 


Statanents  3. 3.  It.  is  initially  void.  Here  it  is  taken  to  be 

the  1 X 1 or  2 X 2 block  starting  at  Lll. 

Statements  3.3.1e.  This  is  the  search  for  alxior2x2  block 
described  in  Section  2. 

Statanent  3.3.1e.3.  The  subroutine  EXCHNG  [n]  is  used  repeatedly 
to  move  the  block  just  located  to  the  beginning  of  A22*  After  each 
exchange  of  a conplex  block  SPLIT  [ii]  is  called  to  recompute  its  eigen- 
values and  to  see  if,  owing  to  rounding  error,  it  can  be  split  into  a pair 
of  real  eigenvalues. 

Statement  3.3.2.  Because  A^^  is  void,  A^^  is  effectively  deflated. 

Statement  3.3.3.  The  matrix  A^2  saved  below  the  lower  subdia- 
gonal of  A,  in  case  the  atten^^t  to  deflate  A^^j^  is  unsuccessful.  Since 
the  routine  SHRSLV  [1],  which  computes  the  deflating  transformation, 
requires  Aj^j^  to  be  lower  Hessenberg,  the  subroutine  DAD  is  called  to 
transform  A^^j^  and  SHRSLV  has  been  modified  to  return  with  a signal 

if  SOTie  element  of  the  deflating  transformation  exceeds  the  bound  RMAX. 
Otherwise  the  matrix  P that  determines  the  transformation  overwrites 
Aj^2*  is  once  again  called  to  restore  A^^j  to  its  original  form. 

Statement  3.3.5.  The  submatrix  A22»  which  was  overwritten  by 
SHRSLV,  must  be  restored  before  another  attenpt  to  deflate  is  made. 

Statement  3.4.  Only  if  A22  is  not  void  was  a deflating  transforma- 
tion generated. 

Statement  3.6.  Set  Lll  to  point  to  the  next  submatrix  before  contin- 


uing. 


ooo  oo  noonnonnoooorinooonnnnonoonoooonnonr.onnonnnr.onr-.nnnooooonnnr:  on 


bUdROlJTiNE  RDIAG  ( A . LOA » N» EPSHOR » RMAX » ER * F.  I » TYPP , BS»  X , LDX , F A ID 
INTEGER  LUAf  LDX.  N»  TYPE(N)»  RS(N) 

WEAL  A(LDA»N>»  ER(N)*  EITn)*  X(LDX»N)»  EPSHfR»  PmAx 
LOGICAL  fail 

HOIAG  REDUCES  A MATRIX  A TO  BLOCK  DIAGONAL  FORM  HY  flPST 
REDUCING  IT  TO  0UASI-TRIAN6ULAR  FORM  BY  HQP3  AMP  THEm  HY 
SOLVING  THE  MATRIX  EQUATION  -All*P+P*A?P=Al2  TO  INTrOOUCE  7ER0S 
ABOVE  THE  DIAGONAL.  THE  PARAMETERS  IN  THE  CAI  LING  SEQUENCE 
are  (STARRFP  PARAMETERS  APE  ALTERED  BY  THE  SUBROUTINE): 

i^A  AN  ARRAY  THAT  INITIALLY  CONTAINS  THE  M X M MATRIX 

TO  HE  REDUCED.  ON  RETURN*  A CONTAINS  THE  REDUCFD 
BLOCK  diagonal  MATRIX. 

LI)A  THE  LEADING  DIMENSION  OF  ARRAY  A. 

N THE  ORDER  OF  THE  MATRICES  A AND  X. 

LPSHQR  THE  CONVERGENCE  CRITERION  FOP  SUBROUTINE  hQR3. 

RMax  THE  MAXIMUM  SI7E  ALLOWED  FOR  ANY  FLFMEfjT  oF  THE 

TRANSFORMATIONS. 

A SINGLY  SUBSCRIPTED  REAL  ARRAY  COMTATMImG  THE  REAL 
PARTS  OF  THE  EIGENVALUES. 

KK I A SINGLY  SUBSCRIPTED  REAL  ARRAY  CONTAING  THE  IMAGINARY 

PARTS  OF  THE  EIGENVALUES. 

►TYpt  A SINGLY  SUBSCRIPTED  INTEGER  ARRAY  WHOSE  I-TH  ENTRY  TS 

0 IF  THE  I-TH  EIGENVALUE  IS  REAL 

1 IF  THE  I-TH  EIGENVALUE  IS  COMPLEX  WITH  POSITIVE 
IMAGINARY  PART 

2 IF  THE  I-TH  EIGENVALUE  IS  COMPLEX  wITh  NEGATIVE 
IMAGINARY  PART 

-1  IF  THE  I-TH  eigenvalue  COULD  NOT  HE  COMPUTED 

*HS  a singly  subscripted  integer  ARRAY  THAT  CONTAINS  BLOCK 

STRUCTURE  INFORMATION.  IF  THERE  IS  A BLOCK  OF  ORDER 
K STARTING  AT  A(L*L)  IN  THE  OUTPUT  MATRIX  A,  THEN 
85(L)  CONTAINS  THE  POSITIVE  INTEGfR  K.  RS(L+l)  CONTAINS 
-(K-l)f  BS(LA+2)  = -(K-2).  ...»  RS(L+K-1)  = -1. 

THUS  A POSITIVE  INTEGER  IN  THE  L-TH  ENTRY  OF  RS 
INDICATES  A NEW  BLOCK  OF  ORDER  BS(L)  STARTING  AT  A(L,D. 

*X  AN  ARRAY  INTO  WHICH  THE  REDUCING  TRANSFORMATIONS 

ARE  TO  BE  multiplied. 

LDx  THE  LEADING  DIMENSION  OF  ARRAY  X. 

►fail  A LOGICAL  VARIABLE  WHICH  IS  FALSE  ON  NORMAL  RETURN  AND 

TRUE  IF  THERE  IS  ANY  ERROR  IN  BDIaG. 


bDIAG  USES  SUBROUTINES  ORTHES.  ORTRAN.  HOR3»  EXCHMG*  SHRSLV.  PAD.  AND  SPl  IT.  ? 

I 

INTEGER  DAll.  DA22.  I*  J*  K*  KMl.  KM2*  L.  LEAVE.  LOOP. 

1 Lll.  L22.  L22M1.  NK  i 

RLAL  C.  CAV.  L).  El.  E2.  RAV.  SC.  TEMP 

FAIL  = .TRUE.  I 

convert  a to  upper  hessenherg  FORM,  I 

call  ORTHES  (LOA.N. l.N. A,ER)  i 
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C 

E 

C 

C 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


call  ORTRAN  (LOAfN»l»N>A»ER*X) 

CONVERT  A TO  QUASI-UPPER  TRIANfiULAR  FORM  BY  OR  mfThOO. 
call  HQR3  (A»X»N.l»N»EPSHQR»ER»El»TYPE»LnA.LDX) 

CHECK  TO  SEE  IF  HQR3  FAILED  IN  COMPUTING  ANY  EIGENVALUE 

DU  5 I = 1»N 

IF  (TYPE(I).EQ.-l)  60  TO  QOO 

continue 


REDUCE  A TO  BLOCK  DIAGONAL  FORM 


SEGMENT  A INTO  4 MftTRICE^: 


All»  A DA 11  X DA  11  block 
A(L11»L11)»  


DA22 


whose  (Irll-ELEMENT  IS  AT  A(Lll»Lll)»  A2?»  A DA??  X 
BLOCK  WHOSE  (1»1)-ELFMFNT  IS  AT  A(L22»L2?)J  Al?, 

A DAll  X DA22  BLOCK  WHOSE  ( 1 » 1 ) -ELEMENT  IS  AT  A(Lll»L??)> 
AND  A21»  A g^22  X DAll  BLOCK  = 0 WHOSE  (Itl)- 


elemenT  is 


:L22>L11 ) . 


This  LOOP  USES  lii  as  loop  index  and  splits  oee  a block 
starting  at  A(L11»L11) . 

Lll  = 1 

10  ASSIGN  550  TO  LOOP 

ASSIGN  600  TO  LEAVE 
IF  (Lll  .GT.  N)  GO  TO  LtAVE 
L22  = Lll 

THIS  LOOP  USES  DAll  AS  LOOP  VARIABLE  AND  ATTEMPTS  TO  SPLIT 
OFF  A BLOCK  OF  SIZE  DAll  STARTING  AT  A(L11»L11) 

lOu  ASSIGN  350  TO  LOOP 

IF  (L22  .NE.  Lll)  GO  TO  lin 
UAll  = TYPE(Lll)  + 1 
L22  = Lll  ♦ DAll 
L22M1  = L22  - 1 
GO  TO  290 
llO  CONTINUE 

COMPUTE  THE  AVERAGE  OF  ThE  EIGENVALUES  IN  All 

RAV  = 0. 

CAV  = 0. 

DO  120  I = L11»L22M1 
HAV  = RAV  + ER(I) 

CAV  = CAV  + ARSCEI ( I ) ) 

120  CONTINUE 

RAV  = RAV  / FLOATiDAll) 

CAV  = CAV  / FlOAT(DA11 ) 

LOOP  ON  eigenvalues  OF  A?2  TO  FIND  THE  ONE  CLOSEST  TO  THE  AVERAGE, 

0 = (RAV-ER(L22) )»*2  + ( CAV-EI ( L22 ) ) **? 

K = L22 

L = L22  + TYPEtL22)  + 1 
130  ASSIGN  145  TO  LOOP 

IF  (L.GT.N)  GO  TO  150 
C = (RAV-ER(L) )**2  + (CAV-EI (L> ) **2 
IF  (C.GE.D)  GO  TO  140 
K = L 


140 

CONTINUE 

L = L ♦ TYPE(L)  ♦ 1 

145 

GO  TO  130 

r 'll 

K ■ 

C 

c 

1 50 

CONTINIJE 

C 

LOOP  TO  MOVE  THE  EIGENVALUE  JUST  l.OCATED 

C 

INTO  FIRST  POSITION  OF  BLOCK  A22. 
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assign  280  TO  LOOP 

IF  (TYPE(K) .NE.O)  GO  TO  200 


THE  BLOCK  WE*RE  MOVING  TO  AHD  TO  All  IS  A 1 X 1 

NK  = 1 
CONTINUE 

IF  (K  .EG.  L22)  GO  TO  280 
KMI  ^ K ^ 1 

IF  (TYPE{KM1 ) .EQ.O)  GO  TO  IQO 

WE'RE  SWAPPING  THE  CLOSEST  BLOCK  WITH  A 2 X 2 
KM2  = K - 2 

CALL  EXCHNG  (A»X»N»KM2.2»1»EPSH0R.FAIL.I  0A,LDX) 

IF  (FAIL)  GO  TO  900 

TWY  TO  SPLIT  THIS  BLOCK  INTO  2 PFAL  EIGENVALUES 

CALL  SPLIT  (A*X»N»KMl»Fl»E2»LnA»LDX) 

IF  (A(K»KM1).EQ.0.)  GO  TO  170 

BLOCK  IS  STILL  COMPLEX. 

TYPE(KM2)  = 0 
TYPE(KMl)  = 1 
TYPE(K)  = 2 
ER(KM2)  = EP(K) 

EI(KM2)  = 0. 

EH(K)  = El 
ERtKMl ) = El 
EI(KMl)  = E2 
EI(K)  = -E2 
GO  TO  180 

COMPLEX  BLOCK  SPLIT  INTO  TWO  REAL  EIGENVALUES. 

CONTINUE 
TYPE(KMl)  = 0 
TYPE (KM2)  = 0 
ER(KM2)  = ER(K) 

ER(KMl)  = El 
ER(K)  = E2 
EI(KM2)  = 0. 

El  (KMl ) = 0. 

K = KM2 

IF  (K  .LE.  L22)  GO  TO  280 
GO  TO  160 

WE'RE  SWAPPING  THE  CLOSEST  BLOCK  WITH  A 1 X 1. 
CONTINUE 

call  EXCHNG  (A»X»N»KM1, 1»1»FPSH0R»FAIL»I  DA,LDX) 

IF  (FAIL)  GO  TO  900 
TEMP  = ER(K) 

ER(K  ) = ER(KMl) 

ER(KMl)  = TEMP 
K = KMl 

IF  (K  .LE.  L22)  go  TO  280 
GO  TO  160 

THE  BLOCK  WE'RE  MOVING  TO  AOO  TO  All  IS  A 2 X 2. 

CONTINUE 
NK  = 2 
CONTINUE 

IF  (K  .EO.  L22)  GO  TO  280 
KMl  = K - 1 

IF  (TYPE(KMl)  .EG.  0)  GO  TO  240 

WE'RE  SWAPPING  THE  CLOSEST  BLOCK  WITH  A 2 X 2 BLOCK, 
KMj2  = K - 2 , 

CALL  EXCHNG  ( A » X» N» KM2 , 2 » 2» EPSHOR » FA IL» LOA , LOX ) 


non  non  non  n non  non  non  non 
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IF  (FAIL)  GO  TO  900 


TRY  TO  SPLIT  SWAPPED  BLOCK  INTO  TWO  PEAI  S. 

CALL  SPLIT  (A.X»N»K»El.E2*LnA»LDX) 

ER(KM2)  = ERJK) 

ER(KMl)  z ER(K+1) 

EKKM2)  z EI(K) 

EI(KMl)  z EKK  + l) 

IF  (A(K+1»K)  .EG.  0.)  GO  TO  220 

STILL  COMPLEX  BLOCK. 

ER(K)  z El 
ER(K+1)  z El 
EI(K)  z £2 
EKK  + 1)  z -E2 
GO  TO  230 

TWO  REAL  ROOTS. 

CONTINUE 
TYPE(K)  z 0 
TYPE(K+1)  z 0 
ER(K)  z El 
ER|K+1)  = E2 
eiIk)  z 0. 

EI(K+1)  z 0. 

CONTINUE 
K = KM2 

IF  (K.EQ.L22)  GO  TO  260 
60  TO  210 

WE»PE  SWAPPING  THE  CLOSEST  BLOCK  WITH  A 1 X 1. 

CONTINUE  , „ 

CALL  EXCHNG  ( A . X » N» KMl , 1 » 2 » EPSHQR f FA IL . I DA , LDy ) 

IF  (FAIL)  GO  TO  900 

TYPE  (KMl)  z 1 

TYPE(K)  z 2 

TYPE(K+1)  z 0 

ER(K  + D z ER(KMl) 

FR(KMl)  z ER(K) 

EI(KMI)  = EI(K) 

EI(K)  = EI(K+1) 

EKK  + 1)  z 0. 

60  TO  2S0 

CONTINUE 
K = KMl 

IF  (K.EO.L22)  GO  TO  260 
GO  TO  210 

TRY  TO  SPLIT  RELOCATED  COMPLEX  BLOCK. 

CONTINUE 

CALL  split  (A»X»N»K»El»E2.LnA»Lnx) 

IF  ( A(K+IfK) .EQ.O. ) 60  TO  270 

STILL  COMPLEX. 

ER(K)  z El 
ER(K+1)  z El 
EI(K)  z E2 
EI(K+1)  = -E2 
GO  TO  280 

SPLIT  into  two  real  EIGENVALUES. 

CONTINUE 
TYPE(K)  z 0 
TYPE(K+1)  z 0 
ER<K)  z El 
FR(K+1)  z E2 
El (K)  = 0. 

EI(K+1)  = 0. 


non  oonr-.o  nnri  on  r.onnno  nnn  nnn  nnnnn  nnn  non 
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?B0 

?Q() 


CONTINUE 

OAll  = DAll  + NK 
L22  = Lll  + DAll 
L22M1  = L22  - 1 
CONTINUE 

ASSIGN  400  TO  LEAVE 
IF  (L22  ,6T.  N)  GO  TO  LEAVF 


attempt  to  SPLIT  OFF  A RLOCK  OF  SIZE  DAll. 


DA22  = N - L22  1 


SAVE  Al2  IN  ITS  TRANSPOSE  FORM  IN  BLOCK  A21. 

no  300  J r L11»L22M1 
DO  300  I = L22»N 
A(I»J)  = A(J.I) 

300  CONTINUE 


COrjVERT  All  TO  LOkVER  QUASI~TRIANGULAR  AND  MULTIPLY  IT  BY  -1  AND  CHAriBF 
A12  appropriately  (FOR  SOLVING  -A1 1*P+P* A2?=ft 12 ) . 

CALL  DAO  ( A.LDA*Lll.L22Ml.Lll»N» 1. *0) 
call  dap  (A»LDA.L11»L22M1.L11»L22M1»-1.,1) 

SOLVE  -All*P  + P*A22  = A12, 

CALL  SHRSLV  ( A ( Ll 1 » L 1 1 ) » A ( L22 » L22 ) . A ( L 1 1 . L22 ) » 

1 UAll  »nA22»LUA.L')A.LDA.RMAX.FAIL) 

ASSIGN  400  TO  LEAVE 
IF  (.NOT.  FAIL)  GO  TO  LEAVF 

CHANGE  All  BACK  TO  UPPER  OUASI-TRI ANGULAR. 

CALL  DAD  (A»LDA.Lll.L22Ml.Lll»L22Ml.l..l  ) 
call  dad  (A»LDA.LllfL22Ml.Lll»L22Ml.-l. ,0) 

WAS  UNABLE  TO  SOLVE  FOR  P - TRY  AGAIN 


MOVE  SAVED  A12  BACK  INTO  ITS  CORRECT  POSITION. 

DO  310  J=L11*L22M1 
DO  310  I = L22»N 
A(J»I)  = A(I*J) 

A(1»J)  = 0. 

310  CONTINUE 


3biJ  GO  TO  100 
4U(i  CONTINUE 

change  SOLUTION  TO  P TO  PROPER  FORM. 

IF  (L22  .GT.  N)  GO  TO  440 

call  dad  ( A.LDA.L11»L22M1»L11»N. 1. .0) 
call  DAO  (A.L0A»L11»L22M1#L11»L22M1»-1.,1) 

MULTIPLY  TRANSFORMATION  INTO  X. 

only  columns  L22  THRU  N ARE  AFFECTED. 

DO  410  U = L22.N 
DO  410  I = 1»N 

DO  410  K = l11»L22M1 

X(I.J)  = X(I»J)  ♦ X(l»K)  * A(K,J) 

4lfi  CONTINUE 

ZERO  OUT  A12  FOR  EASE  IN  HANDLING. 

DO  420  J = L22»N 

DO  420  I = L11*L22M1 

A( If J)  = U. 

420  CONTINUE 


nor-,  ooo  onnn  oo  o 
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430 

44u 

450 

4b0 

470 

480 

500 

510 

52u 

550 

f)Un 

900 


ZEHO  OUT  TRIANGULAR  BLOCK  RELOW  DlAGONAl . 

DO  430  J = Llii^22Ml 
DO  430  I = L22»N 
A(1>J)  = 0. 

CONTINUE 

SCALE  THOSE  COLUMNS  OF  X THAT  WON»T  BE  ALTFREO  AGAIN  TO  Ur'lTY. 
CHANGE  All  APPROPRIATELY. 

CONTINUE 

DO  500  J r L11»L22M1 


88 


450  I = 1»N 
SC  = SC  + (X(I* J) )**2 
CONTINUE 
SC  = SQRT(SC) 

DO  460  I = 1»N 

X( I»J)  = X(I»J)  / SC 
CONTINUE 

DU  470  I = L11»L22M1 
A(I»J)  = A(I.J)  / SC 
CONTINUE 

DO  480  1 = L11»L22M1 
A(J*1)  = A(J,I)  * SC 

continue 

CONTINUE 

STORE  BLOCK  SIZE  IN  ARRAY  BS. 

bS(L11)  = DAll 
J = DAll  - 1 
IF  (J  ,EQ.  0)  GO  TO  520 
DO  510  I = IfJ 

BS(L11+I)  = -(DAll-I) 
CONTINUE 
CONTINUE 
Lll  = L22 
GO  TO  10 
CONTINUE 
FAIL  = .FALSE. 

RETURN 

ERROR  return. 

continue 

fail  = .TRUE. 

return 

END 


- 27  - 


bUdROUTlNE  DAU  (A»  NA,  H»  12,  Jl,  J2,  R,ISW) 

IF  ISW  = 0,  SUBROUTINE  DAn 

COMPUTES  0*A  WHERE  D IS  THE  MATRIX  WITH  ONPS  DOWN 
THE  MINOR  DIAGONAL  AND  A IS  THE  INPUT  MATRIX. 

IF  ISW  = 1,  IT  COMPUTES  THE  PRODUCT  A*0. 

IT  COMPUTES  THIS  PRODUCT  FOR  ROWS  II  THRU  I?  AND  COL- 
UMNS Jl  THRU  J2.  IT  ALSO  MULTIPLIES  EACH  El  EMENT  OF 
THE  PRODUCT  WITH  THE  CONSTANT  R.  NA  IS  THE  FIRST 
DIMENSION  OF  THE  MATRIX  A,  THE  PRODUCT 
OVERWRITES  the  SPECIFIED  LOCATIONS  OF  MATRIX  A. 


REAL  A(NA» 1) , R 

IE  18  E8  !§8 

NRD2  = IFIX( (12  - II  + 1)  / 2 
DO  100  J = J1,JP 
DO  50  IPl  = 1,NRD2 


I = IPl  - 1 
TEMP  = A(I1+I,J) 

A(I1+1,J)  = A(I2-I,J)  ♦ R 
A( I2-I » J)  = temp  * R 
50  CONTINUE 

lou  continue 

IF  (M0D(I2-I1,2)  .EQ.  1)  RETUI 
I = II  + NKD2 

no  CONTINUE 
RETURN 
150  continue 

DO  160  J = Jl,J2 

A(I1,J)  = A(I1,J)  * R 

160  continue 

RETURN 


1)  RETURN 


computes  the  product  ad  where  D is  as  ABOVE. 


continue 

IF  (Jl  ,i 


IF  (Jl  ,E(i»  J2)  60  TO  350 
NCU2  = IFIX( ( J2  - Jl  + 1) 
DO  300  JPl  = 1,NCD2 
J = JPl  - 1 
UO  250  I = I1»I2 
TEMP  = A(I,J1+J) 
A(I,J1+J)  = A(I,J2-J) 
A(I,J2-J)  = TEMP  * R 
25(1  CONTINUE 
SOU  CONTINUE 

IF  ( M0D(J2-J1,2)  .EQ.  1) 
J = Jl  ♦ NCD2 
DO  310  1=11,12 

A(I»J)  = A(I,J)  ♦ R 

310  continue 

RETURN 

350  CONTINUE 

UO  360  I = 11,12 

A(I,J1)  = A(I,J1)  ♦ R 

360  Continue 

RETURN 

END 


return 


-r-'* 
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subroutine  SHRSLV  U»R*C»M.N»NA»NB»MC»RMaX»FAIL) 

SHHSLV  is  a fortran  IV  subroutine  to  solve  the  real  matrix 

EQUATION  AX  ♦ XB  = C»  WH  RE  A IS  IN  LOWER  REAL  SCHljR  FORM 
ANU  0 IS  IN  UPPER  REAL  SCHUR  FORM,  SHRSlV  USES  THE  AUX- 
ILLIARY  SUBROUTINE  SYSSLV»  WHICH  IT  COMMUNICATES  wITh  THROUGH 
THE  COMMON  BLOCK  SLVBLK.  THE  PARAMETERS  IM  THE  CALLING 
SEQUENCE  ARE 

A A DOUBLY  SUBSCRIPTED  ARRAY  CONTAING  THE  MATRIX  A IN 

LOWER  SCHUR  FORM 

H A DOUBLY  subscripted  ARRAY  CONTAINING  THE  MATRIX  R 

IN  UPPER  REAL  SCHUR  FORM 

C A DOUBLY  subscripted  ARRAY  CONTAINING  THE  MATRIX  C. 

M THE  ORDER  OF  THE  MATRIX  A 

N THE  ORDER  OF  THE  MATRIX  B 

NA  THE  FIRST  DIMENSION  OF  THF  ARRAY  A 

NH  THE  FIRST  DIMENSION  OF  THE  ARRAY  F3 

IlC  THE  FIRST  DIMENSION  OF  THE  ARRAY  C 

NMAX  MAXIMUM  ALLOWED  SIZE  OF  ANY  ELEMENT  OF  THE  TRaNSFORMATTON 
KAIL  INDICATES  IF  SHRSLV  FAILED 

INTEGER  M»N»NA,NB»NC»K»KMl.riK»KK»L»LMl  ,DL»Ll  » I » IR»  J*  JA.NSYS 
RLAL  A(NA>l)f  M(NBrl)»  C ( NC » 1 ) • T»  P 
logical  FAIL»  SING 

COmMON/SLVBLK/  T(5»5) »P(b) ,NSYS»SING 
FAIL  = .TRUE. 

L = 1 

) LMl  = L - 1 
DL  = 1 

IF  (L  .EO.  N)  GO  TO  IS 
IF  (B(L+1»L)  .NE.  0.)  DL  = P 
j LL  = L+DL-1 

IF  (L  ,EQ.  1)  GO  TO  3U 
DO  20  J = L»LL 
DO  20  I = 1»M 


pu  continue 

30  K = 1 

40  KMl  = K - 1 
DK  = 1 

IF  (K  .EQ.  M)  GO  TO  4b 
IF  (A(K»K+1)  .NE.  0.)  DK  = P 

4b  KK  = K+DK-1 

IF  (K  .EO.  1)  GO  TO  60 
DO  bO  I=K»KK 
DO  50  J=L»LL 
00  50  JA=1.KM1 

C(I*J)  = C(I»J)  - A(I»JA)*C(JA» J) 

50  continue 

60  IF  (DL  .EO.  2)  GO  TO  80 
IF  (DK  ,FQ.  2)  GO  TO  70 
T(l»l)  = A(K»K)  ♦ B(L»L) 

IF  (T(l»l)  .EQ.  0.)  RETURN 
C(K»U  = C(K»L)  / T(l.l) 

IF  (AbS(C(K»D)  .GE.  RMAX)  RETURN 
GO  TO  100 

70  T(l»l)  = A(K.K)  ♦ B(L.L) 

T(l»2)  = A(K»KK) 

T(2*l)  = A(KK.K) 

T(2»2)  = A(KK»KK)  B(L»L) 

P(I)  = C(K»L> 


DO  20  IP  = 1*LM1 
C(I»J)  = C(I.J) 


- C( I . IB)*B( IP» J) 


I tu. 


P(2)  = C(KK»L) 

NSYS  = 2 
CALL  SYSSLV 
IF  (SING)  RETURN 
C(K»L)  = P(l) 

IF  (AHS(C(KK»L) ) 
GO  TO  100 


RMAX)  RETURN 
, RMAX)  RETURN 


IF  (DK 
T(l»l) 
T(l»2) 
T(2»l) 

irm 

P(2)  = 
NSYS  = 

call  S 
IE  Tsi 


DK  .EQ*2)  GO 

1)  = A(K»K)  * 

2)  = H(LL»L) 

1)  = B(L»LL) 

2)  = A(KtK)  ♦ 
= C(K»L) 

= C(KfLL) 

= 2 

SYSSLV 

SING)  RETURN 


TO  90 
B(LfL) 


B(LL»LL) 


C(K»L)  = P(l) 

IF  (ABS(C(K»D) 
C(K»LL>  = P(2) 

IF  (AHS(C(K*LL) ) 
GO  TO  100 


RMAX)  RETURN 
, RMAX)  RETURN 


T (1»1) 
TU»2) 
T(l»3) 
T(l»4) 
T(2»l) 
T(2»2) 
T(2*3) 
r(2»4) 
T(3»l) 

T(3»4) 
T(4»1) 
T(4»2) 
T14»3) 
T (4»4) 
P(l)  = 
P12)  = 
P(3)  = 
P(4)  = 
NSYS  = 


= A(K»K)  + 
= A(K.KK) 

= B(LL*L) 

= 0. 

= A(KK»K) 

= A(KK»KK) 
= U. 

= T(l»3) 

= H(L»LL) 

= 0. 

= A(K»K)  ♦ 
= T(l»2) 

= 0. 

= T(3»l) 

= T(2.1)  , 
= A(KK»KK) 
C(K»L) 
C(KK.L) 

C (K»LL) 


B(L»L) 


B(L.L) 


B(LL»LL) 


B(LL.LL) 


P(4)  = C(KK*LL) 
NSYS  = 4 
call  SYSSLV 
IF  (SING)  RETURN 
C(K»L)  = P(l) 

IF  (ABS(C(KfD)  . 
C(KK,L)  = P(?) 

IF  ( ABS(C(KK*L) ) 
C(K»LL)  = P(3) 

IF  (ABS(C(K»LL) ) 
C(KKtLL)  = P(4) 

IF  ( ABS(C(KK.LL) ) 
K = K + DK 
IF  (K  .LE.  M)  GO 
L = L + DL 
IF  (L  .LE.  N)  GO 
FAIL  = .FALSE. 
return 

END 


RMAX)  RETURN 
. RMAX)  RETURN 
, RMAX)  RETURN 


RMAX)  return 


TO  40 


i 
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SUHROUTINE  SYSSLV 

bYsbLV  IS  A FORTRAN  IV  SUBROUTINE  THAT  SOLVES  THE  LINEAR 
system  AX  = B of  ORDER  N LESS  THAN  5 BY  CROUT  REDUCTION 
FOLLOWED  0Y  BLOCK  SUBSTITUTION.  THE  MATRIX  A.  THE  VECTOR  H. 

AND  THE  ORDER  M ARE  CONTAINED  IN  THE  ARRAYS  A.R»  AtjD  THE 
VARIABLE  N OF  THE  COMMON  BLOCK  SlVBLK,  THE  SOLUTION  IS  RETURmED 
IN  THE  ARRAY  B. 

common  /SLVBLK/  A(5*5)»  B(5)»  N»  SING 
REAL  MAX 

logical  sing 

SING  = .TRUE. 

I NMl  = N - 1 
N1  = N + 1 

compute  The  lu  factorization  of  a 

DO  80  K=1»N 
KMl  = K-1 

IF  (K  .EO.  1)  GO  TO  20 
L)0  10  I=K»N 
DO  10  J=1»KM1 

A(I.K)  = A(IfK)  - A(I.J)*A(J»K) 

10  CONTINUE 

?0  IF  (K  .EG.  N)  GO  TO  100 

KPl  = K+1 
MAX  = ABS( A(K»K) ) 

INTR  = K 
DO  30  I=KPl»N 

AA  = AdS(A( I»K) ) 

IF  (AA  ,LE.  MAX)  GO  TO  30 
MAX  = AA 
INTR  = I 
30  CONTINUE 

IF  (MAX  .EQ.  0*)  RETURN 
A(N1»K)  = INTR 
IF  (INTR  .EO.  K)  GO  TO  SO 
DO  40  J=1.N 
TEMP  = A(K»J) 

A(K»J)  = A(INTR.J) 

A(INTR»J)  = TEMP 
40  CONTINUE 

SO  DO  bO  J=KP1»N 

IF  (K  .EQ.  1)  60  TO  70 
DO  80  1=1. KMl 

A(K.J)  = A(K.J)  - A(K.I)*A(I.JI 

60  CONTINUE 

70  A(K.J)  = A(K.J)  / A(K.K) 

BO  CONTINUE 

INTERCHANGE  ThE  COMPONENTS  OF  B 


GO  TO  110 


C 

C 

C 


lOO  DO  110  Jrl.NMl 

INTR  = A(Nl.J) 

IF  (INTR  .E(3.  J) 

TEMP  = B(J) 

B(J)  = B(INTR) 

B(INTR)  = TEMP 
no  CONTINUE 

SOLVE  SY  = B 

?U0  B(l)  = B(l)  / A(l.l) 

DO  220  1=2. N 
IMl  = I-l 
DO  210  J=1.IM1 

B(I)  = B(I)  - A(I.J)*B(J) 
210  CONTINUE 

H(I)  = B(I)  / A(I»I) 

220  CONTINUE 

SOLVE  UX=Y 


■■ill  Hi 


300  UO  310  II=1»NM1 
I 

L)h  3i5^J=I1»N 

8(1)  = B(I) 
310  continue 

blNG  = .FALSE. 

KtTURN 

END 
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c 

c 

c 

c 

c 
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c 

c 

c 

c 

t 

c 

c 

c 

c 
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c 
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SUHHOUT INE  H0R3  ( A » V . N . NLOW , MIJP , EPS  * ER  p F 1 , T YPE  » N A , Nv  ) 

INTEGER  NpNApNLOWpNUPpNVp  TYPE<N) 

HtAL  A(NApN)  pEKN)  pER(N)  pEPSpV(NVpN) 

H(JH3  REOKCES  THE  UPPER  HESSENBFR6  MATRIX  A TO  OiiaSi- 
TKIANGULAH  FORM  BY  UNITARY  SIMILARITY  TRANSFORMATIONS, 
TmE  EIGENVALUES  OF  Ap  m/HICH  ARE  CONTAINED  IN  THr  1X1 
ANU  2X2  DIAGONAL  BLOCKS  OF  THE  REDUCED  MaTRIXp  aPE 
OKUEREl)  in  descending  order  of  magnitude  ALONG  THE 
diagonal.  the  transformations  are  ACCUMULATED  IN  THE 
ARRAY  V.  HOR3  REQUIRES  ThE  SUBROUTINES  FXCHNGp 
QRSTFPp  and  SPLIT.  ThE  PARAMETERS  IN  THr  CALLING 
sequence  are  (STARRED  PARAMETERS  ARF  ALTERED  HY  THE 
SUBROUTINE) 


♦ A 


*V 


N 

NLOW 

NUP 


EPS 

*ER 

*EI 

• TYPE 


NA 

NV 


AN  ARRAY  THAT  INITIALLY  CONTAINS  THE  N X M 
UPPER  HESSENBLRG  MATRIX  TO  RE  REDUCED.  ON 
RETURN  A CONTAINS  THE  RFDUCEOp  QUASI- 
triangular  matrix. 

AN  ARRAY  THAT  CONTAINS  A MATRIX  INTO  WHICH 
THE  REDUCING  TRANSFORMATIONS  ARE  TO  PE 
MULTIPLIED. 

THE  order  of  the  MATRICES  A AND  V. 

a(nlow-Ipnlow)  and  a(nuppnup+u)  are 

ASSUMED  To  RE  ZERO.  AND  ONLY  ROWS  N|  OW 
THROUGH  NUP  AND  COLUMNS  NLOW  THROUGH 
NUP  ARE  TRANSFOPMEDp  RESULTING  IN  FhE 
CALCULATION  OF  EIGENVALUES  NLOW 
THROUGH  NUP. 

A CONVERGENCE  CRITERION. 

AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE  REAL 
PARTS  OF  THE  EIGENVALUES. 

AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE 
IMAGINARY  PARTS  OF  THE  EIGENVALUES. 

ANU  INTEGER  ARRAY  WHOSE  I-TH  ENTRY  TS 

0 IF  the  I-TH  eigenvalue  is  RFAl  p 

1 IF  THE  I-TH  eigenvalue  IS  COMPLEX 
with  POSITIVE  IMAGINARY  PART. 

2 IE  THE  i-TH  eigenvalue  IS  COMPLEX 
with  negative  imaginary  PART. 

-1  IF  THE  1-th  eigenvalue  WAS  NOT 
CALCULATED  SUCCESSFULLY. 

the  " ‘ 

THE 


first  dimension  of  thf  array  a. 
first  dimension  of  the  array  V. 


c 

c 

c 

c 

c 

c 

c 


internal  variables 

INTEGER  Ip  ITpLpMUpNL»NU 
REAL  E 1 p E2  p P p p R p S p T » w p X p y p ? 

LOGICAL  FAIL 

INITIALIZE. 

DO  10  irNLOWpNUP 
TYPE(I)  = -1 
If)  CONTINUE 
T = 0. 

MAIN  LOOP.  FIND  AND  ORDER  EIGENVALUES. 

NU  = NUP 

lUO  IF(NU  .LT.  NLOW)  60  TO  500 
IT  = 0 

OR  LOOP.  find  NEGLIGABLE  ELEMENTS  AND  PEPFORm 
UR  STEPS. 

110  CONTINUE 

SEARCH  BACK  FOR  NEGLIGARLE  ELEMENTS. 

L = NU 

120  CONTINUE 


onn  ooo  noon  non  oooo  nono 


- 33  - 


IF(L  .EG.  NLOW)  GO  TO  I3n 

IF(ABS(A(L»L-1) ) .LT.  EPS* ( ABS ( A (L-1 »L-1 ) ) +ABS ( A ( L»L> ) ) > 
GO  TO  130 
L = L-1 
GO  TO  120 
CONTINUE 

TEST  TO  SEE  IF  AN  EIGENVALUE  OR  A ?yi?  BLOCK 
HAS  BEEN  FOUND, 

X = A(NU,NU) 

IF(L  ,EQ.  NU)  GO  TO  300 
Y = A(NU-1*NU-1) 


TEST  iteration 
IT  IS  10  OR  20 


OUNT,  IF  IT  IS  30  QUIT, 
ET  up  an  AD-HOC  SHIFT, 


IF(IT  ,EQ,  60)  GO  TO  500 
IF(IT,NE,40  ,AND,  1T,NE,50)  GO  TO  150 

AO-HOC  SHIFT, 

T = T ♦ X 
UO  140  I=NLOW»NU 

Ad.I)  = A(I*I)  - X 
CONTINUE 

S = AHS(A(NU,NU-1) ) + AHS ( A ( NU-1 , NU-2 ) ) 

X = 0,75*S 

Y = X 

M = -n,4375*S**2 

CONTINUE 
IT  = IT  ♦ I 

LOOK  FOR  T*0  CONSECUTIVE  SMALL  SUB-DIAGONAL 
elements, 

NL=  NU-? 

CONTINUE 

Z = A(NL»NL) 

K = X - Z 
S = Y - Z 

P = (p*S-W)/A(NL+l*NL)  + A(NL»NL+1) 

Q r A(NL+1»NL+1)  - Z - R - S 
K = A(NL+2»NL+1) 

S = AHS(P)  ♦ ABS(Q)  ♦ APSIR) 

P = P/S 
G = Q/S 
H = R/S 

IF(NL  ,EG,  L)  60  TO  170 

IF(ABS(ATNL*NL-l)i*(APS(Q)+ABS(P) ) ,LF, 
EPS*ArtS (P) * ( AHS ( A ( ML-1 f NL-1 ) ) +ABS ( 7 ) +Af 
GO  TO  170 
NL  = NL-1 
60  TO  IbO 
CONTINUE 

PERFORM  A QR  STFP  BETWEEN  NL  AND  NU. 

CALL  QRSTEP ( A » V r P , 0 » R » NL  » NU » N » N A » NV ) 

GO  TO  110 

2X2  BLOCK  FOUND. 

IF(NU  ,NE.  NLOW+1)  A(NU-1»NU-2)  = 0. 

A(NU,NU)  = A(NU»NU)  ♦ T 
A(NU-1#NU-1)  = A(NU-1»NU-1)  ♦ T 
TYPE(NU)  = 0 
TYPE (NU-1)  = 0 
MU  = NU 

LOOP  TO  POSITION  2X2  BLOCK, 

CONTINUE 
NL  = MU- I 


APS<A(NL4l,NL+l ) ) ) ) 


non  non  nnnn  non  non 
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C ALL  E XCHNG ( A » V » N » MU » 1 » 2 » EPS » F A I L » NA  * M V ) 
IF(.NOT.  FAIL)  GO  TO  32b 
TYHE(MU)  = -1 
TYP 
TYP 
GO  TO  500 
CONTINUE 
MU  = MU+2 
GO  TO  340 

330  CONTINUE 

THE  NEXT  BLOCK  IS  1X1. 

IF(AHS(A(MU*MU) ) .GE.  ARS ( A ( MU+ 1 » MU+ J > ) ) 
1 GO  TO  350 

CALL  EXCHNG( A»V»N»MU» 1 , 1 . EPS » F A IL » NA * MV ) 
MU  = MU+1 

340  CONTINUE 

GO  TO  320 
35(1  CONTINUE 

MU  - NL 
NL  = 0 

IF(MJ  .NE.  0)  GO  TO  310 

GO  hack  AND  GET  THE  NEXT  EIGENVALUE. 

4U0  CONTiNUf: 

NU  = L-1 
GO  Tn  100 


(MU+1)  = -1 
(MU+2)  = -1 


ALL  THE  FIGNVALUES  HAV/E  BEEN  FOUND  ANn  ORDERFfi. 
COMPUTE  THEIR  VALUES  AND  TYPE. 


500  IP(NU  .LT.  NLOW)  GO  TO  507 
DO  5U3  1=1 » MU 

A(I»I)  = A(I»I)  ♦ T 
503  CONTINUE 
507  continue 
NO  - NUP 

510  continue 

IF(TYPE(NU)  .ME.  -1)  GO  TO  515 
NU  = NU-1 
GO  TO  540 
515  CONTINUE 

IF(MU  .EG.  NLOW)  GO  TO  520 
IF(A(NU»NU-1 ) .EQ.  0.)  GO  TO  520 


520 


0.)  GO  to  520 


2X2  HLOCK. 

CALL  SPLIT (A»V»M»NU-1. El »E2»NA.NV) 
IF(A(NU»NU-1)  .EO. 

ER(NU)  = El 
El  (NU-1)  = E2 
ER(NU-l)  = ER(NU) 

El  (NU)  = -EI(NU-l) 

TYPE(NU-l)  = 1 
TYPE(NU)  = 2 
NU  = NU-2 
GO  TO  530 
CONTINUE 

SINGLE  ROOT. 


ER(NU)  = A(NU»NU) 
EKNU)  = 0. 

NU  = NU-1 
530  CONTINUE 
540  CONTINUE 

IF(NU  .GE.  NLOW)  GO  TO  510 

RETURN 

END 


non  nnn  non  non  r.oo  nnor-,  noon  noo  o 
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1 

2 

3 


215 

22(1 


1 

2 


22b 

230 

300 


310 


32o 


E 

C 


attempt  to  split  the  block  into  two  real 
eigenvalues. 

CALL  SPLIT(A. V*MfNL»El»E2»NA»NV) 

IF  THE  SPLIT  WAS  SUCCESSFUL#  GO  AND  ORDER  ThF 
REAL  EIGENVALUES. 

IF(A(MU»MU-1)  .EQ.  0.)  GO  TO  310 

TEST  TO  SEE  IF  THE  BLOCK  IS  PROPERLY  POSITIONEDt 
ANO  IF  NOT  EXCHANGE  IT 

IF(MU  .EG.  NUP)  GO  TO  400 
IF(MU  .EQ.  NUP-1)  GO  TO  220 
IF(A(MU+2»MU+1)  .EO.  0.)  60  TO  220 

THE  NEXT  BLOCK  IS  2X2. 

IF( A(MU-1,MU-1 )*A(MU»MU)-A(MU“1»MU)*A(MU#MU-1) 

. 6E  . A ( MU+ 1 » MU+ 1 ) * A ( Ml  1+2 » MU+  2 ) - A ( MU+ 1 f Ml  )+2 ) * 
A(MU+2#MU+1) ) 

GO  TO  400 

CALL  E XCHNG (A»V#N#NL#?»2» EPS # FA I L » MA , MV ) 

IFl.NOT.  FAIL)  GO  TO  215 
TYPE(ML)  = -1 
TYPE(NL+1)  = -1 
TYPECML-^2)  = -1 
TYPF(NL+3)  = -1 
GO  TO  500 
CONTiNUf 
MU  = MU+2 
GO  TO  230 
CONTINUE 

THE  NEXT  BLOCK  IS  1X1. 

IF  < A ( MH-1  » MU- 1 ) * A ( MU  » MU ) -A  ( MU- 1 » MU  ) * A ( MU  » Ml)- 1 ) 

,6E.  A(MU+1»MU+1)**2) 

GO  TO  400 

CALL  EXCHNGIA#  V»N»NL»?» 1 .EPS# FAIL# NA. MV) 

IF( .NOT.  FAIL)  GO  TO  225 
TYPE(ML)  = -1 
TYPE(NL-Kl)  = -1 
TYPE(NL+2)  = -1 
GO  TO  500 
CONTINUE 
MU  = MU+1 
CONTINUE 
GO  TO  210 

SINGLE  EIGENVALUE  FOUND. 

NL  = 0 

A(NU#NU)  = A(NU#NUJ  + T 

IF(NU  .NE.  NLOW)  A(NU#NU-1)  = 0. 

TYPE(NU)  = 0 
MU  = NU 

LOOP  TO  POSITION  ONE  OR  TWO  REAL  EIGENVALUES. 

CONTINUE 

POSITION  THE  EIGENVALUE  LOCATED  AT  A(NL#NL). 

CONTINUE 

IFIMU  .EQ.  NUP)  GO  TO  350 
IFIMU  .EQ.  NUP-1)  GO  TO  330 
IF(A(MU+2»MU+1)  .EQ.  0.)  GO  TO  330 

THE  NEXT  BLOCK  IS  2X2. 

IF(A(MU#MU)**2  .GE. 

A ( MU+ 1 # Ml  )•♦■  1 ) *A  ( MU+2  # MU+2 ) - A ( Mil+ 1 , MU+? ) ♦ A ( MU+2 » Ml  l ) ) 
GO  TO  400 


I IL.  . 
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SUhROUTINE  EXCHMGC A*V»N»L»H1»R?»EPS»FAIL.NA»NV) 

INTEGER  Hl»B2.L»NAtNV 
real  A(NA»N) *EPS»V(MV»N) 

logical  fail 

given  the  upper  HESSENBERG  r^ATRlX  A WITH  COMSECIiTiVE 
HlXHl  AND  B2XB2  DIAGONAL  BLOCK*;  (Hl»R?  .LE.  2) 

STARTING  AT  A(LtL)f  EXCHNG  PRODUCES  A UNITARY 
SIMILARITY  transformation  THAT  EXCHANGES  THF  BLOCKS 
ALONG  WITH  THEIR  EIGENVALUES.  THE  TRANSFORMATION 
IS  ACCUMULATED  IN  V.  EXCHNG  HfOUIRES  THF  SURROllTiNE 
ORSTEP.  the  PARAMETERS  IN  THE  CALLING  SEQUENCE  ARE 
(STARRED  parameters  ARE  ALTERED  BY  THE  SUBROUTINE) 

*A  THE  MATRIX  WHOSE  BLOCKS  ARE  TO  PE 

INTERCHANGED, 

*V  THF  ARRAY  INTO  WHICH  THE  TRANSFORMATIONS 

ARE  TO  RE  accumulated. 

N THE  ORDER  OF  THE  MATRIX  A. 

L THE  POSITION  OF  THE  BLOCKS. 

Fll  THE  SIZE  OF  THE  FIRST  BLOCK. 

H2  THE  SIZE  OF  THE  SECOND  BLOCK. 

EPS  A CONVERGENCE  CRITERION, 

♦FAIL  A LOGICAL  VARIABLE  WHICH  IS  FALSE  ON  A 

NORMAL  RETURN.  IF  THIRTY  ITERATIONS  WpRE 
performed  without  CONVERGENCE*  FAIL  IS  SET 
TO  TRUE  ANU  the  ELEMENT 
A(L+B2.L+B2-1)  CANNOT  BF  ASSUMED  ZERO. 

NA  THE  FIRST  DIMENSION  OF  THE  ARRAY  A. 

NV  THE  first  DIMENSION  OF  THE  ARRAY  V. 

INTERNAL  VARIABLES, 

INTEGER  I*IT,J»L1*M 
REAL  P,Q,R*S,W,X, Y,Z 

FAIL  = .FALSE. 

IF(H1  ,E(i.  2)  GO  TO  4i) 

IF(B2  .EG.  2)  GO  TO  10 

INTERCHANGE  1X1  ANU  1X1  BLOCKS. 

LI  = L+1 

0 = A(L+1.L+1)  - A(L.L) 

P = A(L»L+1) 

R = AMAXKP.Q) 

IF(R  .EG.  0.)  return 
P = P/R 
0 = 0/K 

H = SORT(P^^2  ♦ Q**2) 

P = P/R 
0 = 0/K 
DO  3 J=L»N 

S = P^A(LtJ)  ♦ Q*A(L+1»J) 

A(L+1»J)  = P^A(L+1»J)  - G^A(L.J) 

A(L.J)  = S 
3 CONTINUE 

DO  5 I=1*L1 

S = P^Ad.L)  + Q*A(I,L+1) 

A(I*L+1)  = P^A(I»L+1)  - O^Ad.L) 

Ad*L)  = S 
b CONTINUE 

DO  7 1=1, N 

S = P^V(  I ,L)  + Q*Vd  ,L  + 1 ) 

Vd,L  + l)  = P^V(I,L  + 1)  - Q*Vd,L) 

Vd,L)  = S 
7 CONTINUE 

A(L+1»L)  = 0. 

RETURN 
10  CONTINUE 

INTERCHANGE  IXI  AnD  2X2  BLOCKS. 


I 


■x±. 


ooo 


i 


X = A(L»L) 

P = 1. 

Q = 1. 

P = 1. 

CALL  QHSTEP (A»V»P»Q»H»L» L+2 * N » NA » MV ) 

^0  IT  = ?T+1 

IF(IT  ,LE.  60)  GO  TO  30 
FAIL  = .TRUE. 

RFTURN 
30  CONTINUE 

P = A(L»L)  - X 
0 = A(L+1»L) 

^ aCl'^Srstep  (A»V.P*Q.R,L.  L+2  » N , NA  , NV  ) 

IF(ABS(A(L+2»1 +1) ) .GT. 

1 FPS»(AHS(A(L  + 1»L-»-1)  )+ABS(A(L+2»L  + ?)  ) ) ) 

I GO  TO  20 

A(L+2»L+1)  = 0. 

RETURN 
40  CONTINUE 

INTERCHANGE  2X2  AND  B2XB2  BLOCKS. 

H = L+2 

IF(H2  .eg.  2)  M = M+1 
X = A(L+1»L+1) 

Y = A(L»L) 

W = A(L+1.L)»A(L»L+1) 

P = 1. 

Q = 1. 

H = 1. 

call  QRSTeP( A» V»P»0»R»L»M.N»NA*NV) 

IT  = 0 

50  IT  = IT+1  ^ ^ 

IF(IT  .LE.  60)  GO  TO  60 
FAIL  = .TRUE. 
return 
on  CONTINUE 

7 - A(LfL) 

R = X - 2 
5 = Y " Z 

P = (H*S-W)/A(L+1»L)  + A(L.L+1) 

Q = A(L+1.L+1)  - Z - H - S 
K = A(L+2fL+l)  , . 

S = ABS(P)  + ABS(Q)  + ABS(R) 

P = P/S 
0 = 0/S 
K z P/S 

CALL  QRSTEP(A»V»P»Q»R*L»M.N»NA»NV) 

IFTah5(A(m-1,M-2) ) .GT.  EPS*(AbS(A(M-l ,M-1 ) )+AH5(A(m-2»M-?) ) ) ) 

1 GO  TO  50 

A(M-1*M-2)  = 0. 

return 

continue 

END 
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subroutine  split ( A» V »N»L. El »E2,NA»NV) 

INTEGER  L»N»NA*NV 
REAL  A(NA»N) f V(NV»N) 


given  the  upper  HESSENbERb  MATRIX  A WITH  A oLOcK 

starting  at  a(l»l)»  split  determines  if  the 
corresponding  eigenvalues  are  real  or  complex,  if  they 
ARE  REAL*  A ROTATION  IS  DETERMINED  THAT  REDUCES  THE 
BLOCK  TO  yPP  ‘ - 

OF  LARGEST  A _ 

ROTATION  IS  accumulated  IN  V 
OR  COMPLEX)  ARE  RETURNED  IN  El 
IN  THE  CALLING  SEQUENCE  ARE  TS 
ALTERED  BY  THE  SUBROUTINE) 


TRIANGULAR  FORM  WITH  Thf  FlGFNVAl  UE 
OLUTE  value  appearing  FIhST.  THE 


THE  EIGENVALUFS  (HEaL 
AND  E2.  THE  PARAMETERS 
TARRED  PARAMETERS  ARE 


*A  THE  UPPER  HESSEMVERG  MATRIX  WHOSE  2X? 

BLOCK  IS  TO  BE  SPLIT. 

♦V  The  ARRAY  IN  WHICH  THE  SPLITTING  TRANS- 

FORMATION IS  TO  BE  ACCUMULATED. 

N The  order  of  the  matrix  a. 

L The  position  of  thf  2X2  BLOCK. 

♦El  ON  RETURN  IF  THE  EIGENVALUES  ARF  COMPLEX 

♦E2  El  CONTAINS  THEIR  COMMON  REAL  PART  AND 

E2  contains  The  positive  imaginary  part. 

IF  THE  EIGENVALUES  ARE  REAL.  El  CONTAINS 

The  one  largest  in  absolute  value  and  F2 

CONTAINS  THE  OTHER  ONE. 

NA  The  first  dimension  of  thf  array  a. 

NV  the  first  dimension  of  the  array  V. 


internal  VARIABLES 


Integer  i.u.li 

REAL  P.O.R.T.U.W.X.Y.Z 


X = A{L+1»L+1) 

Y = A(L.L) 

w = A(L.L+1)TA(L+1.D 
P = (Y-X)/2. 

0 - P**2  + W 

IF(Q  .GE.  0.)  GO  TO  b 

complex  eigenvalue. 

El  = P ♦ X 
E2  = SORT{-Q) 

RETURN 

5 continue 

Two  REAL  EIGENVALUES.  SET  UP  TRANSFORMATION. 
L - SORT(O) 

IF(P  .LT.  0.)  GO  TO  10 
Z = P + Z 
GO  TO  20 
10  CONTINUE 

Z,=  P - z 
20  continue 

IF(Z  .EO.  0.)  GO  TO  30 
R = -W/Z 
GO  TO  40 
30  CONTINUE 
R = 0. 

40  continue 

IF(ABS(X+Z)  .GF.  ABS(X+R))  7 = R 

Y = Y - X - Z 
X = -Z 

T = A(L.L+1) 

U = A(L+1.L) 

IF(ABS(Y)+ABS(U)  .LE.  ABS ( T ) +AhS ( X ) ) GO  TO  GO 

^ = V 

GO  TO  70 


I 


faO  continue 

Q = X 
P = T 

70  continue 

K = SQKT(P**2  + Q**2> 
IF(R  .GT.  0. ) GO  TO  80 
El  = A(L»U 
E2  = A(L+1*L+1) 
A(L+1»L)  = 0. 

RETURN 
80  CONTINUE 
P = P/R 
0 = Q/R 

premultiply. 


DO  90  J=L»N 
2 = A(L*J) 

A(L»J)  = P*Z  + Q*A(L+1»J) 
A(L+1*J)  = P*A(L+1*J)  - Q*Z 
CONTINUE 

POSTMULTIPLY. 

LI  = L+1 
DO  100  I=1*L1 
2 = A(I»L) 

A(I*L)  = P*Z  + 0«A(l*L+l) 
A(I»L+1)  = P*A(I»L+U  “ Q*2 

continue 

ACCUMULATE  THE  TRANSFORMATION 

DO  110  1=1 »N 
2 = V(I»L) 

V/(  I »L)  = P*2  + Q*V(  I »L+1 ) 
V(I*L+1)  = P*V(I»L+l)  - (i*Z 

continue 

A(L+1»L)  = 0. 

El  = A(L»L) 

E2  = A(L+1»L+1) 

RETURN 

END 
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SUHHOUT INE  0»STEP ( A » V r P » 0 » R , NL , NU»  N»  NA , MV ) 

INTEGER  N»NA»NLf NU»MV 
REAL  ^ (NA»N) fP»Q»R>V(NVf N) 

OHSTFP  performs  one  implicit  or  step  ON  THF 
UPPER  HESSEN0ER6  MATRIX  A.  THF  SHIFT  IS  DETFRMTMtD 
HY  THE  NUMbERS  P*Q»  AND  K.  ANU  THE  STFP  IS  APPLIED  TO 
ROwS  AND  COLUMNS  NL  THROUGH  NU.  THE  TRA^JSFOPMATIONS 
ARE  ACCUMULATED  IN  V.  THE  PARAMETERS  IN  THF  CA|  lING 
sequence  are  (STARRED  APRAMFTERS  ARE  ALTFREO  PY  THf 
SUBROUTINE) 


The  upper  hessenberg  matrix  on  which  thf 

OR  STEP  IS  TO  HE  PERFORMED. 

The  array  in  which  the  transformations 
ARE  TO  PE  ACCUMULATED 

parameters  that  determine  the  SHIFT. 


The  LOWER  limit  of  the  step. 

THE  UPPER  LIMIT  OF  THE  STEP. 

THF  ORDER  OF  THE  MATRIX  A. 

THE  FIRST  dimension  OF  THE  ARRAY  A. 
the  first  dimension  of  THE  ARRAY  V. 


INTERNAL  VARIABLES. 

INTEGER  I»U»K.NL2»NL3»NUMI 
REAL  S»X.Y»Z 
LOGICAL  LAST 

NL2  = NL+2 
DO  10  I=NL2»NU 
AiI»I-2)  = 0* 
in  continue 

IF(NL2  .EQ.  NU)  GO  TO  30 
NL3  = NL+3 
UO  20  I=NL3»NU 
A(I»I-3)  = 0. 

20  CONTINUE 
30  continue 
NUMl  = NU-1 
UO  130  K=NLfNUMl 

determine  THE  TRANSFORMATION. 

LAST  = K .EO.  NUMl 
IF(K  .EQ.  NL)  GO  TO  40 
P = A(K»K-1) 

Q = A(K+1»K-1) 

R = 0. 

IF(. NOT. LAST)  R = A(K+2»K-1) 


X = ABS(P)  ♦ ABS(O)  + APS(R) 

1F(X  .EQ.  0.)  GO  TO  130 
P = P/X 
Q = 0/X 
R = H/X 
CONTINUE 

S = SQRT(P**2  ♦ Q**2  + R**2) 

IF(P  ,LT.  0.)  S = -S 
IF(K  .FO.  NL)  60  TO  50 
A(K»K-1)  = -S*X 
60  TO  faO 
CONTINUE 

IF(NL  .NE.  1)  A(K.K-l)  = -A(K»K-1) 
CONTINUE 
P r P + S 
X = P/S 
Y = Q/S 
Z = R/S 
a - Q/P 
R = R/P 


♦ Q**f 

S = 

60  TO 


+ APS(R) 
130 


R**2) 


PREMULTIPLY 
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DO  8i)  J=K»N  . ,, 

P = a(K>J)  + Q*A(K+1»J) 

IF  (LAST)  60  TO  7() 

p = P R*A(K  + 2»  J)  , _ , 

A(K+2*J)  = A(K+?»J)  - P*Z 

+ = A(K+1.J)  - P*Y 

A(K.J)  = A(K»J)  - P*X 
CONTINUE 

POSTMULTIPLY. 

J = MIN0(K+3»NU) 

^^P^z^xiAn^K)  + Y*A(I.K  + 1> 
IMLAST.^SOJ0^90^ 

A(I»K+2)  = A(I»K+2)  - P*R 

A?lTK+n  = A(I»K+U  - P*Q 
A(I»K)  = A(I»K)  - P 
CONTINUE 

accumulate  the  transformation  in 


UO  120  I=1»N 

P = X*\i  ( I »K  ) + 
IF(LAST)  GO  TO 


Y*V(IrK+l) 

110 


130  continue 
RETURN 
FNU 


- p.R 

viiIk+u  = VlI.Ktll  - P*o 
V( I»K)  = V(I»K)  - p 
CONTINUE 


■c, 


J 
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subroutine  0RTHES(NM»N»L0W» 1GH,A»0RT) 

1 NTEGER  I»J»MiN*lI>UJ*LArMP»NM,IGH»KPl»LOM 
REAL  A(NM»N) »0RTCIGH) 


+ ARS(A( ) 
0.)  GO  TO  IflO 


REAL  A(NM»N) »0RTCIGH) 

REAL  F»6*H»SCALE 
LA  = IGH  - 1 
KPl  = LOW  1 
1F(LA  .LT.  KPl)  GO  TO  200 
DO  180  M=KP1»LA 
H = (J. 

ORT(M)  = 0. 

SCALE  = 0, 

UO  90  1=1 » IGH 
SCALE  = SCALE  + ARS(A(I»M 
IF(SCALE  .EQ.  0.)  GO  TO  1 
MP  = M + IGH 
DO  100  II=M»IGH 
I = MP  - II 

ORT(I)  = A(ItM-l)/SCALE 
H = H ORT(I)*ORT(I) 
CONTINUE 

G = -SIGN(SQRT(H) .OHT(M) ) 
H = H-  ORT(M)*G 
ORT(M)  = ORT(M)  - G 
DO  130  J=M»N 
F = 0. 

DO  110  II=M»IGH 
I = MP  - II 
F = F + ORT(I)*A(I»J) 
CONTINUE 

IF  (H  .EO.O. ) GO  TO  lb2 
F = F/H 

DO  120  I=M.I6H  . 

A(I»J)  = A(I»J)  - F*0 
CONTINUE 

continue 

DO  160  I=1»IGH 
F = 0. 


.J)  - F*0RT(I) 


DO  140  JJ=M»IGH 
J = MP  - JJ 
F = F ORT(J)*A(I» J) 

140  CONTINUE 

F = F/H 

DO  150  J=M»IGH 

A(I»J)  = A(I.J)  - F*ORT(J) 
150  CONTINUE 

160  CONTINUE 
162  CONTINUE 

ORT(M)  = SCALE*ORT(M) 

A(M»M-1)  = SCALE*G 
180  CONTINUE 

200  return 

END 


- 43  - 


oO 

BU 


SUBROUT I ME  ORTR AN ( NM » N » LOW » 1 6H , A » ORT ; ? ) 
INTEGER  I » J»N»KL*MM»MP»NM» IGH*L0W»MP1 
HLAL  A(NM» IGH) *ORT(IGH) »Z(NM»N) 

REAL  G*H 
no  80  I=1»N 
UO  60  J=1>N 
Z(1»J)  = 0. 


CONTINUE 
= 


Z<I»I)  = 1 

continue 

KL  = IGH  - LOW  - 1 ^ ^ 


MP  = IGH  - MM  , 

H r a(MP»MP-1)*ORT(MP) 
IF(H  .EQ.  0.)  60  TO  140 

1 Ol) 

VIPI  = MP+1 

DO  100  I=MP1»I6H  . 

ORT(l)  = A(I»MP-1) 
CONTINUE 

UO 

UO  130  J=MP»I6H 
G — 0 • 

DO  110  1=MP.IGH 

a z 6 ORT(I)*ZU*J) 

CONTINUE 

1 Zo 

G - G/H 

DO  1?0  I=MP»I6H 

Z(I»J)  = Z(I»J)  + 6*0RT 
CONTINUE 

130 

CONTINUE 

140 

continue 

?un 

return 

END 
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ABSTRACT  fContlnua  on  rovoroo  olda  H nacaooary  and  IdantUy  by  block  numbar) 


This  paper  describes  an  algorithm  for  reducing  a real  matrix  A to  block 
diagonal  form  by  a real  sLidlarity  transfoimation.  "Hie  columns  of  the  trans 
formation  corresponding  to  a block  span  a reducing  subspace  of  A,  and  the 
block  is  the  representation  of  A in  that  subspace  with  respect  to  the 
basis.  The  algorithm  attenpts  to  control  the  condition  of  the  transforma- 
tion matrices,  so  that  the  reducing  subspaces  are  well  conditioned  and  the 
basis  vectors  are  numerically  independent. W 


edition  of  I NOV  6S  IS  OBSOLETE 


UNCLASSIFIED > 

security  classification  of  this  page  rwhen  neim  Enlrrei 


